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Simulating protostellar jets simultaneously at launching and 

observational scales 

Jon P. Ramsey and David A. Clarke^ 
ABSTRACT 

We present the first 2.5-D MHD simulations of protostellar jets that include 
both the region in which the jet is launched magnetocentrifugally at scale lengths 
< 0.1 AU, and where the propagating jet is observed at scale lengths > 10'^ AU. 
These simulations, performed with the new AMR-MHD code AZEuS, reveal inter- 
esting relationships between conditions at the disc surface, such as the magnetic 
field strength, and direct observables such as proper motion, jet rotation, jet 
radius, and mass flux. By comparing these quantities with observed values, we 
present direct numerical evidence that the magnetocentrifugal launching mecha- 
nism is capable, by itself, of launching realistic protostellar jets. 

Subject headings: magnetohydrodynamics — ISM: jets and outflows — stars: 
formation 



1. Introduction 

Jets and outflows from protostellar objects are fundamental aspects of the current star 
formation paradigm, and are observed anywhere star formation is ongoing. The mechanism 
proposed by Blandford & Payne (1982), in which jets are launched from accretion discs 
by gravitational, magnetic, and centrifugal forces, has been extensively studied numerically 
{e.g., Uchida & Shibata 1985; Meier et al. 1997; Ouyed & Pudritz 1997a,b, 1999; Krasnopol- 
sky et al. 1999; Vitorino et al. 2002; von Rekowski et al. 2003; Ouyed, Clarke, & Pudritz 
2003; Porth & Fendt 2010; Staff et al. 2010). By treating the accretion disc as a boundary 
condition {e.g., Ustyugova et al. 1995), one can study jet dynamics independently of the disc 
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{e.g., Pudritz et al. 2007) though, in order to resolve the launching mechanism, numerical 
simulations have not followed the jet beyond 100 AU {e.g., Anderson et al. 2005). 

In stark contrast, protostellar jets are > 10^ AU long (Bally, Reipurth, & Davis, 2007), 
and only recently have observations reached within 100 AU of the source {e.g., Hartigan, Ed- 
wards, & Pierson, 2004; Coffey et al, 2008). This large scale difference between observations 
and simulations makes direct comparisons difficult and, in this work, we aim to close this 
gap. We present axisymmetric (2.5-D) simulations of protostellar jets launched from the in- 
ner AU of a Keplerian disc, and follow the jet well into the observational domain (2500 AU). 
These calculations allows us to address the efficacy of the magnctocentrifugal mechanism, 
and to relate conditions near the disc with directly observable properties of the jet. 

The simulations presented herein are performed with an adaptive mesh refinement 
(AMR) version of ZEUS-3D (Clarke, 1996, 2010) called AZEuS (Adaptive Zone Eulerian 
Scheme). The ZEUS-3D family of codes are among the best tested, documented, and most 
widely used astrophysical MHD codes available, though this is the first attempt to couple 
ZEUS-3D with AMR^ We have implemented the block-based method of AMR detailed in 
Berger & Colella (1989) and Bell et al. (1994). Significant effort was spent minimising errors 
caused by passing waves across grid boundaries, which is of particular importance to this 
work. A full description of the code and the changes required for AMR on a fully-staggered 
mesh will appear in Ramsey & Clarke (in preparation). 

2. Initialisation 

Observationally, the inner radius of a protostellar accretion disc, ri, is between 3-5 
(Calvet et al., 2000) and, for atypical T Tauri star (M = 0.5 Mq, = 2.5Rq), n = 0.05 AU. 
Thus, following Ouyed & Pudritz (1997a), we initialise a hydrostatic, force-free atmosphere 
surrounding a 0.5 Mq protostar coupled to a rotating disc with = 0.05 AU. However, unlike 
Ouyed & Pudritz we use an adiabatic equation of state that conserves energy across shocks 
rather than an isentropic polytropic equation of state, as the distinction becomes important 
for supermagnetosonic flow (Ouyed, Clarke, & Pudritz, 2003). 

We solve the equations of ideal MHD^ (7 — 5/3) over a total domain of 4096 AU x 



^ENZO, a hybrid N-body Eulerian code (O'Shea et al., 2004), links AMR with the hydrodynamical portion 
of ZEUS-2D. 

^AZEuS solves either the total or internal energy equation. We chose the latter because positive-definite 
pressures trump strict conservation of energy in these simulations; see Clarke (2010). 
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256 AU. To span the desired length scales, nine nested, static grids (refinement ratio 2) are 
initialised each with an aspect ratio of 4:1 (16:1 for the coarsest grid only) and bottom left 
corner at the origin. Our finest grid has a domain 4 AU x 1 AU and a resolution Az = ri/S = 
0.00625 AU which we find sufiicient to resolve the launching mechanism. Thus, the effective 
resolution for the entire domain is > 26 billion zones. The simulation highlighted in §3 was 
run to t = 100 yr with an average time step in the finest grid of ~ 3 minutes and thus ~ 18 
million time steps. 

During the simulations, a thin region of low velocity and high poloidal magnetic field, 
-Bp = y^Bf+B^., develops along the symmetry axis, the edge of which is defined by a large 
gradient in the toroidal magnetic field, drB^. Insufficient resolution of drB^ can lead to 
numerical instabilities, and grids are added dynamically whenever this gradient is resolved 
by fewer than five zones. 



2.1. The atmosphere 

The atmosphere is initialised in hydrostatic equilibrium (HSE; Vz = Vr = = 0). 
Because the LHS of the equation governing HSE, 

Vp + pV0 = O, (1) 

is not a perfect gradient, differencing it directly on a staggered-mesh can commit sufficient 
truncation error to render the atmosphere numerically unstable. Thus, we replace V0 with 
the corresponding poloidal gravitational acceleration vector, 

^=--Vph, (2) 
Ph 

where ph and ph are the hydrostatic density and pressure given by: 

Here, p; and are the initial density and pressure at and p (x p"' \s assumed throughout 
the atmosphere at t = 0. In this way, differencing equation (1) maintains HSE to within 
machine round-off error indefinitely. 

However, equations (3) as given arc singular at the origin where truncation errors are 
significant regardless of resolution. These errors can launch a supersonic, narrow jet from 
the origin destroying the integrity of the simulation. To overcome this problem, we replace 
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the point mass at the origin with a uniform sphere of the same mass and a radius Rq, thus 
modifying the first of equations (3) to: 



X 7-1 
Ph ' 

pi 



r' + z'>Ro'; 



^j,2 _|_ ^2 ' 

(4) 

I p oW^ ' + <Ro ■ 

If Rq is sufficiently resolved {e.g., four zones), the numerical jet is eliminated. The resulting 
"smoothed potential" is superior to a "softened potential" since the former has no measurable 
effects beyond Rq. Here, we use Rq = r^. 

The atmosphere is initialised with the force-free magnetic field used by Ouyed & Pudritz 
(1997a): 

Bi y/r^ + {z + z^y -{z + z^) 



A 



-Dz = , -Dr = 7. — , ^<p — U, 

r or az 



where A^ is the vector potential, is the disc thickness (set to ri), and Bi is the magnetic 
field strength at ri, given by: 

Here, pi and f3i (plasma beta at ri) are free parameters. 

Finally, to ensure the dechning density and magnetic field profiles do not fall below 
observational hmits, we add fioor values pfloor ~ 10~^pi and -B^^floor ~ 10~^Bi (c.f. Bergin & 
Tafalla 2007, Vallee 2003) to equations (4) and (5). By imposing HSE and the adiabatic gas 
law at f = 0, a fioor value on p imposes effective fioor values on g and p as well. 



2.2. Boundary Conditions 

In the accretion disc (-z < 0, r > ri), v^p — vk — \^GM^/r, the Keplerian speed, and 
— C'^K = 10~^vk is an "evaporation speed" at the disc surface. The disc and atmosphere 

— * 

are initially in pressure balance with a density contrast rj — Pdisc/Patm = 100, while B is 
initialised using equations (5). 

Following Krasnopolsky et al. (1999), p, p, and Vz are held constant, Vr — VzBr/B^, 
Vip — vk + VzB^/Bz, Ez{—z) — Ez{z) (where E — v x B is the induced electric field), 
Er{0) = vkBz{0), Er{-z) = Er{0) - Er{z), E^{0) = 0, and E^{-z) = -E^{z). Since is 
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sub-slow, these conditions are formally over-determined and p should probably be allowed to 
float. Indeed, we allowed p to be determined self- consistently in test simulations, and found 
only minor quantitative differences in the jet since the pressure gradient is only about 1% of 
the net Lorentz force at the disc surface. However, allowing p to float in the boundary caused 
undue high temperatures in the disc, and thus small time steps. Therefore, the simulation 
proceeds more rapidly but otherwise virtually unchanged when p is maintained at its initial 
value. 

Inside rj (-2 < 0), we apply reflecting, conducting boundary conditions ( J = V x 7^ 0). 
Thus, p, p, and v are reflected across 2; = 0, and magnetic boundary conditions arc set 
according to Ez{—z) = —Ez{z), Er{—z) = Er{z), and E^{—z) = E^{z). At z = 0, E^ and 
E^ are evolved using the full MHD equations. 

Finally, we use reflecting boundary conditions along the r = symmetry axis, and 
outflow conditions along the outermost r and z boundaries. 



where R is the spherical polar radius. Prom equations (6), (7), and the ideal gas law {p — 
pkT/{m), where (m) is half a proton mass), we derive the following scaling relations to 
convert from unitless to physical quantities: 



2.3. Scaling Relations 
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for 7 = 5/3. Note that is the only free parameter varied in this work. 
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75 km/s (top panel) 
50 km/s (otherwise) 



Fig. 1. — Nested images of a /3i = 40 jet at t = 100 yr. Colours indicate temperature, white 
contours magnetic field lines, maroon contours the slow surface, and arrows the velocity. 
Dashed lines denote grid boundaries. 

3. Results for A = 40 

Figure 1 depicts a jet with /3i = 40 at t ~ 100 yr from the highest resolution grid near 
the disc surface (bottom panel) to the coarsest grid in which the jet has reached a length of 
just under 2500 AU (top panel) ^. A few features worth noting include: 

• When 6 < 60° (angle between Bp and disc surface), Blandford & Payne (1982) show 
that cold gas near the disc is launched into a coUimated outflow. Here, 6 < 60° 
for all r > n, but significant outflow is limited to inside the point where the slow 
surface intersects the disc (rj j ~ 30 AU = jet radius at the disc; second panel from 
top). Below Tj d, cold disc material has moved onto the grid and accelerated into the 



•^Time-lapse animations are available at http://www.ica.smu.ca/zeus3d/rcl0/. 
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outflow. Above rj^d, the weak magnetic field has yet to drive enough disc material onto 
the grid to displace the hot atmosphere, and outflow is stifled. While rj j gradually 
increases with time, the majority of mass flux originating from the disc is driven within 
Ti < r < lOri (0.5 AU; bottom panel). 

• Jet material becomes super-fast (Mf < 5) within a few AU of the disc, and the bound- 
ary between jet and entrained ambient material is defined by a steep temperature 
gradient (contact discontinuity; second panel). Portions of the original atmosphere, 
which remain virtually stationary throughout the simulation, are still visible above 
and ahead of the bow shock (top panel). 

• At large distances from the disc (> 500 AU; top panel), the dynamics of the jet become 
dominated by S^, and the jet is led by an essentially ballistic, magnetic "nose-cone" 
with a Mach number of ~ 10 {e.g.. Clarke, Norman, & Burns, 1986). Still, is a 
small fraction (10~^) of B^, consistent with Hartigan et al. (2007). 

• The knots dominating the bottom panel (c./., Ouyed & Pudritz 1997b) are produced 

— * 

by the nearly harmonic oscillation of Sp in ri < r < 2 ri, whereby 6 fiuctuates between 
55° and 65° with a period ~ 30 Ti. These oscillations result from the interplay between 
in-falling material along the symmetry axis, and under/over pressurisation near the 
central mass. The knots are denser and hotter than their surroundings, and bound by 
magnetic field loops. They occupy a region within ~ 2 AU of the symmetry axis, and 
gradually merge to form a continuous and narrow column of hot, magnetised material^ 
(third panel). As such, they are unhkely to be the origin of the much larger-scale knots 
observed in some jets {e.g., HHlll; Raga et al., 2002). 

Further details of this and other simulations of protostellar jets are left to a future paper, 
and we focus here on a few properties directly comparable with observations. 

4. Comparing simulations and observations 

Table 1 summarises a few observational characteristics of protostellar jets. To connect 
these attributes to conditions in the launching region, we have performed a small parameter 
survey in ^i, and made numerical measurements of the quantities in Table 1. Variation of 
other parameters (such as C, and pi) is left to future work. 



^The knots are resolved by 10-20 zones when they merge, and thus their merger is unlikely related to the 
ever- decreasing resolution of the nested grids. 
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proper motion (km s ^) 100 - 200 (500 max.) 
rotational velocity (km s"^) (5 - 25) ± 5 

FWHM jet width (AU) 30 - 80 (at 200 AU) 
mass-loss rate (10"^ Mq yr'^) 0.01 - 1 

Table 1: Selected observational characteristics of protostellar jets. References: Reipurth & 
Bally (2001), Ray et al. (2007), McKee & Ostriker (2007). 



Note that Pi is the initial value of the plasma beta at n, and not the average (3 in the 
jet. Indeed, Fig. 2a demonstrates that at very early time, {(3) = 8n{p)/{B'^) < (3i/5, where 
B"^ = Bp + B"^, and where the average is taken over zones that exceed a certain threshold 
so that only out-flowing jet material is considered. Thus, the magnetic field within the jet is 
stronger than /3i would suggest. Initially, (p) is dictated by Bp, but becomes dominated by 
B^ within < 10 yr after launch. As time progresses, (P) gradually increases but never rises 
above unity (at least for t < 100 yr), even for » 1. Still, one might speculate from Fig. 
2a that with sufficient time, {/3) — >• 1 regardless of /3i. 



4.1. Proper motion 

For t > 10 yr, the velocity of the tip of the jet, Vjet, is nearly constant and, from Fig. 
2b and Table 2, we find Vj^t oc 

To understand this result physically, we begin with the magnetic forces: 

^11 = -^V|| (rB,) ; 

= ^V|| (rB^) ; (13) 

F± = -^Vx (rB^) + J^Bp, 

{e.g., Ferreira 1997; Zanni et al. 2007) where Vy, V_l are the gradients parallel and perpen- 
dicular to Bp. For a given field line, a stronger Bp at its "footprint" in the disc (r = Tq) 
generates a stronger B^ which leads to stronger gradients in rB^p and thus, from equations 
(13), greater magnetic forces to accelerate the flow. In practice, we flnd that most of the 
acceleration occurs before the fast point (and not the Alfven point) located at r = rf, where 
Tf is a weak function of the field strength at the footprint and thus of B[. 



Following Spruit (1996), one can show that as a function of the "fast moment arm" 
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Fig. 2. — (a) {(5) as a function of time for different (b) Vjet (diamonds) and {v^p) (triangles) 
of each jet as a function of B^. Best fit power-law coefficients for these data are a — 0.44±0.01 
(vjet, solid line) and 0.66 ± 0.01 {{v^), dashed line). 



(^ = Tf/ro), the poloidal velocity at the fast point is: 

since Op^f, the poloidal Alfven speed at the fast point, is roughly proportional to Bi. vk,o = 
^ GM^Itq is the Keplerian speed at the footprint of the field line. We note that measured 
values of Wp^f in our simulations vary as B\-^ and agree with equation (14) to within 1% so 
long as the fluid is in approximate steady-state^. 

After the poloidal force given by equations (13) decreases to 1% of its maximum value (> 
a few Tf), t^p still follows a power law in B\ with index 0.52 ± 0.04 and essentially unchanged 
from equation (14). Nearer the head of the jet where steady state is no longer vahd, we flnd 
(fp) oc ^o-^sio-o^ (-where the momentum- weighted average is taken across the jet radius), only 
shghtly shallower than equation (14). Thus, while the conditions in the jet have changed, 
some memory of the steady-state conditions at rf persists. 

Finally, Vjet (Fig. 2b and Table 2) is within ~ 10% of (wp) near the bow shock and 
maintains the same power- law dependence on B-^. Thus, these jets are essentially ballistic. 



^Indeed, all four steady-state functions from Spruit (1996) remain constant in our simulations to within 
< 5% along steady-state field lines, which we take as validation of our numerical methods. 
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Vjet (km s"^) 


84 


125 


161 


230 


270 


330 


460 


0.44 ±0.01 


{v^) (km s"^) 


2.6 


3.0 


6.2 


10.1 


13.1 


18.4 


31 


0.66 ±0.01 


2rjet (AU) 


21 


40 


60 


85 


94 


104 


130 


0.35 ±0.04 


Mjet (IO-^Mq yr-1) 


0.44 


1.9 


2.8 


4.2 


6.9 


10.1 


17.9 


0.92 ±0.09 



Table 2: Simulation "observables" v^et and (v^) are asymptotic values while rjet and Mjet are 
measured at 2; = 200 AU and t — 20 yr. Uncertainties in a are from the fitting procedure. 



where the observed jet speed fjct oc i^o-^^io.oi ^ short, all measures of jet speed increase 
with Bi, a trend that agrees with Anderson et al. (2005) who find for much less evolved jets, 



4.2. Toroidal velocity 

Figure 2b and Table 2 show averaged over time and the jet volume for z >100 AU as 
a function of B\. Like fjct, v^p asymptotes to a constant value. The region inside 100 AU is 
ignored because the torsion Alfven wave at low z has a non-negligible f;^, is not part of the 
jet, and skews our results. By fitting a power law to these data, we find {v^) oc ij.o-66±o.oi 

Unlike fjct, we have not uncovered a rationale for this power law, yet it seems plausible 
one must exist given the tightness of fit. Eliminating B^ from the power laws for (v^) and 
i^jet, we find that {v^p) oc Uj^f^^^'^^. To render this a useful observational tool, further work is 
needed to quantify the effects of other initial conditions such as C, and pi on both the power 
law index and the proportionality constant, as well as the effect our simplified disc model 
may have on conditions in the jet at observational length scales. 



4.3. Jet radius and mass flux 

The jet radius, rjet, is defined by the contact discontinuity (steep temperature gradient 
in the second panel of Fig. 1) between shocked jet and shocked ambient material, which in 
turn is determined by where the radial jet ram pressure balances all external forces. Since 
ram pressure increases with i>p, rjet should increase with Si, just as observed in Table 2. At 
any given time, we find that rjet varies with B^ as a reasonable power law though, unlike 
i'jet or {vip), the power index is not constant and decreases slowly in time, while rjet itself 
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increases in time, though at an ever-slowing rate. 

The mass flux transported by the jet, Mjet, consists of material from both the disc and 
the atmosphere. Unlike previous simulations where jets are typically evolved long after the 
leading bow shock has left the grid, no part of any bow shock in our simulations reaches 
the boundary of the coarsest grid. Thus, each jet continues to entrain material from the 
atmosphere throughout the simulation at a rate that has a strong dependence on B[, as 
seen in Table 2. Indeed we find that Mjct varies with B[ as a reasonable power law, with 
the power index decreasing slowly in time. As the atmosphere is depleted, the mass fiux 
contribution from the disc (which, by design, is independent of Bi) becomes more important 
and the dependence of Mjet on Bi diminishes. 

5. Discussion 

We have presented the first MHD simulations of protostcUar jets that start from a well- 
resolved launching region (Azmin = 0.00625 AU) and continue well into the observational 
domain (2500 AU). On the AU scale, each jet shows the characteristic and near steady- 
state knotty behaviour first reported by Ouycd & Pudritz (1997b), though the origin of our 
knots is quite different. On the 1000 AU scale, each jet develops into a ballistic, supersonic 
(8 < M < 11) outflow led by a magnetically conflned "nose-cone" (Clarke, Norman, & 
Burns, 1986) and a narrow bow shock, consistent with what is normally observed. 

On comparing Tables 1 and 2, our simulations comfortably contain virtually all observed 
protostellar jets on these four important quantities. We note that these tables would not 
have been in agreement had we stopped the jet at, say, 100 AU and measured these values 
then. It is only because our jets have evolved over five orders of magnitude in length scale 

that we can state with some confidence that the magnetocentrifugal launching mechanism 
is, by itself, capable of producing jets with the observed proper motion, rotational velocity, 
radius, and mass outflow rate. Indeed, our jets are still very young, having evolved to only 
100 yr, and allowing them to evolve over an additional one or two orders of magnitude in 
time may still be useful. For example, it would be interesting to know whether rises 
above unity for any of the jets (Fig. 2a), and thus enter into a hydrodynamically dominated 
regime. It would also be interesting to see how long it takes for the power laws in jet radius 
and mass flux as a function of B-^ to reach their asymptotic limits. 

Our jet widths tend to be higher than those observed, particularly when one considers 
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that the values for rjet in Table 2 are at i = 20 yr^, and that rjet continues to grow in 
time {e.g., for the (3^ = 40 jet, 2 rjet ~ 100 AU hj t = 100 yr). As our jet radii mark the 
locations of the contact discontinuity while observed radii mark hot, emitting regions, our 
widths should be considered upper limits. That our values contain all observed jet widths 
is a success of these simulations. 

Similarly, our numerical mass fluxes are higher than observed values by at least an 
order of magnitude. Since observed mass-loss rates account only for emitting material {e.g., 
in forbidden lines; Hartigan, Morse, & Raymond 1994), and thus temperatures in excess of 
10^ K (Dyson & Williams 1997; p. 104), our mass fluxes arc necessarily upper limits as well. 
Indeed, if we measure our mass fluxes near the jet tip (instead of at 200 AU for Table 2) and 
restrict the integration to fluid above 10^ K, our mass fluxes drop by a factor of 10-100, in 
much better agreement with Table 1. 

We thank the referee for timely and helpful comments on the manuscript, Marsha Berger 
for her AMR subroutines, and Sasha Men'shchikov for early work on AZEuS. Use of MPFIT 
by C. B. Markwardt and JETGET by J. Staff, M. A. S. G. J0rgenson, and R. Ouyed is 
acknowledged. This work is supported by NSERC. Computing resources were provided by 
ACEnct which is funded by CFI, ACOA, and the provinces of Nova Scotia, Newfoundland 
& Labrador, and New Brunswick. 



^Some simulations had not reached t = 100 yr at the time of this writing. 
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